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1.    Introduction 

In  doing  numerical  weather  forecasting  there  are  essentially  two 
problems:   1)  the  forecasting,  which  amounts  to  solving  a  coupled  set  of 
non-linear  partial  differential  equations  in  space  and  time,  and  2) 
determining  the  initial  conditions  for  those  equations  from  observational 
data.   The  problem  of  determining  the  initial  conditions  accurately  is  a 
very  serious  one  which  has  not  yet  been  approached  by  finite  element 
methods.   Thus  we  shall  not  discuss  it  in  this  report. 

The  forecast  equations  are  essentially  the  laws  of  conservations  of 
mass,  momentum  and  energy.   Written  in  the  usual  form,  as  a  set  of  first 
order  partial  differential  equations,  they  are  called  the  "Primitive 
Equations"  or  the  "P.E.  model."  As  such  they  are  hyperbolic,  and  non- 
linear.  These  equations  are 


(tt  +  V'V  +  w  |~)V  +FkxV  +  -Vp  =   vV2V  (la) 
9t                                dZ                                           p 

g  +  if=o  (lb) 

4+^v  +  Miblose  =  ^F  (lc) 

■ff-  +  V-  (pV)  +  f-  (pw)  =  0  (Id) 

dt  dZ 


where  F  is  the  Coriolis  "parameter"   2  ft  sin  9,   0  is  the  potential 

"*/Cp 
temperature   (0  =  T(P/P  )     ) ,   V  is  the  horizontal  velocity  vector 

and  w  is  the  vertical  velocity.   The  wave  solutions  which  result 

consist  of  two  types  of  waves  -  the  so-called  inertial-gravity  waves  and 

the  Rossby  waves  (See  [1]).   It  is  known  that  the  gravity  waves  should 

not  be  present,  and  in  any  numerical  computation  they  die  out,  with  time 
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constants  of  12-24  hours.   However,  this  length  of  time  is  the  time  of 
primary  interest  for  forecast  purposes. 

The  laws  of  conservation  can  also  be  written  in  terms  of  scalar  and 
vector  potentials  i.e.,  a  geopotential  <J>,   a  vorticity  E,     and  a 
stream  function  ij>  .   These  functions  satisfy  a  second  order  elliptic 
and  a  single  hyperbolic  equation  and  thus  have  the  property  that  the 
gravity  waves  are  damped  out,  or  filtered.   (See  equations  (27)  and  (29)). 
This  form  of  the  equations  is  referred  to  as  the  vorticity,  or  filtered, 
form  of  the  equations.   However,  involving  a  second  order  equation  they 
require  different  boundary  conditions,  which  conditions  are  somewhat 
contrived. 

For  the  primitive  equations  (1) ,  one  can  set  up  an  order  of  magnitude 
analysis  and  determine  which  terms  are  significant  for  meteorological 
purposes.   This  leads  to  a  hierarchy  of  equations  which  one  can  solve 
more  easily  than  (1) .   The  zero  order  analysis  gives  the  following  set 

F  k  x  V  =  -  -  Vp  (2a) 

P 

V-(pV)  =  0,  (2b) 

which  dominate  all  the  other  terms  by  an  order  of  magnitude.   This 
suggests  that  the  basic  flow  is  horizontal,  divergence  free,  and  driven 
by  the  balance  of  the  pressure  force  and  the  Coriolis  force. 

At  the  next  level  we  have  the  so  called  barotopic  equations 

(|r  +  V-V)V  +  F  k  x  V  +  V<*>  =  0  (3a) 

a  t 

(|^  +  v-v)<j>  +  <j)V-v  =  o  (3b) 
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where  <j>  =  gh  is  the  geopotential  at  the  height  h,  obtained  by  inte- 
grating —  dp  in  the  z  direction,  since  equations  (3)  are  independent 
of  vertical  variables. 

If  we  take  a  perturbation  approach  about  some  uniform  velocity  V 
and  a  height  H,  we  get  a  set  of  hyperbolic  first  order  equations 

iM=-Vo.Vu+Fv-fi  (4a) 

|^-^.VV-Fu_|£  (4b) 

where  V  -  iu  +  jv  .   If  we  look  for  a  solution  of  the  form 

.    ,  ,v    .   i(wt-k»x) 

Tp  =  (u,v,<j>)  =  i|>  e  (5) 

we   find  three  waves 

w     =  V  •   k  (6a) 

o  o 

W         =  V   •   k  ±VgH  k«k  +  F  (6b) 


The  first  is  a  convective  wave,  due  to  the  uniform  flow.   The  second  and 

third  are  "inertial-gravity"  waves.   If  the  terms   (Fv,  -  -~)     were  zero 

in  (4a)  and   (~Fu,  7~)  were  zero  in  (4b),  then  these  twc  waves  would 

dy 

disappear.   These  are  just  the  terms  in  (2a),  which  essentially  "balance 
out." 

If  the  initial  conditions  contain  any  of  the  eigenvectors  to  which 
these  last  two  waves  are  the  eigenvalues,  then  the  linear  solution 
contains  them.   Thus  any  solution  to  the  non-linear  problem  which 
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initially  contains  some  component  in  these  "directions"  will  propagate 
them  for  some  time.   The  real  problem  is  non-linear,  and  in  actuality 
they  die  out  but  we  would  hope  that  they  are  not  there  to  start  with. 

As  a  result  there  is  considerable  interest  in  getting  correct 
initial  conditions  from  the  observational  data.   We  turn  first  to  the 
methods  which  involve  vorticity.   These  all  involve  the  assumption  that 
there  is  a  static  balance  of  the  wind  and  pressure  fields,  and  that  one 
can  be  derived  from  the  other.   Proceeding  from  there  is  complicated 
by  the  fact  that  in  the  tropics  there  are  very  small  pressure  changes, 
so  that  observational  errors  in  pressure  can  be  as  great  as  the  pressure 
changes  themselves,  while  the  winds  tend  to  be  more  accurately  recorded. 
The  converse  is  true  in  the  mid  -  and  high  -  latitudes.   So  in  the 
mid-latitudes  one  takes  the  pressure  measurements  at  various  points 
(invariably  not  grid  points),  computes  the  pressure  at  the  grid  points 
by  some  weighted  average,  converts  to  the  geopotential  <j>,   and  then 
computes  the  stream  function  ty     using  one  of  the  three  relations  (9)  - 
(11).   The  winds  are  then  obtained  from  \\>.      In  the  tropics  the 
proceedure  is  reversed.   The  observed  winds  V  are  used  to  compute  the 
vorticity 

K   =  k-V  x  V  (7) 

The  stream  function  \\)     is  obtained  from  the  Poisson  equation, 

V2^  =  g  (8) 

Then  the  geopotential  tj>  is  obtained  by  solving  one  of  the  following: 

v2<j>  =  fv2^  (9) 
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V2<fr   =   V(FV^)  (10) 

V2*   =   2J(||  ,  |i)   +  V.(FV^)  (11) 

where     F     is   again   the   Coriolis   factor, 

F  =   20,  sin9  (12) 

J  is  the  symbol  for  the  Jacobian,  and   9   is  the  latitude  from  the 
equator. 

These  equations  are  called  the  quasi-geostrophic  relation,  the 
linear  balance  equation,  and  the  non-linear  balance  equation.   Note  from 
(12)  that  to  use  (9) ,  (10)  or  (11)  to  solve  for  $  at  the  equator 
involves  solving  a  singular  differential  operator.   Note  also  that  these 
equations  are  all  Poisson  equations.   Note  also  that  this  procedure  only 
deals  with  one  part  of  the  wind,  the  rotational  part.   Nowhere  does  the 
divergent  part  enter  the  calculation.   Since  the  wind  is  assumed 
non-divergent  in  the  first  place,  that  is  not  essential  at  the  start. 
However,  a  finite  difference  scheme  will  introduce  some  errors,  and  this 
may  be  significant. 

2.    Finite  Difference  Models 

These  methods  were  used,  and  are  still  used  in  some  places,  for 
many  years.   However,  more  and  more  people  are  turning  to  the  primitive 
equations.   So  let  us  consider  the  simplest  version,  the  barotropic 
model.   First  we  note  that  the  gravity  waves  (6b)  have  a  much  higher 
speed  (about  300  meters/second)  than  the  meteorological  mode  (6a)  (5 
meters/second).   For  computing,  using  a  system  of  hyperbolic  equations, 
the  Courant-Lewy-Fredericks  criterion  requires  that  the  time  step  be 
less  than  Ax  divided  by  the  wave  speed.   These  fast  waves  then  dominate 
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the  size  of  the  time  step,  for  calculating  using  an  explicit  calculation. 
Thus  there  is  considerable  interest  in  semi-implicit  schemes. 

Williamson  and  Browning  [2]  run  some  computer  tests  on  the  primitive 
equations,  using  the  barotropic  model,  writing  the  equations  in  spherical 
coordinates.   They  use  both  the  advection  form,  (a  is  the  Earth's  radius) 


3u       u   3u   v  3u  ,  .-    ,    u  ,    _N  3d     ...  oN 

It-  =  "  IZo^  n  ~  1 1e  +  (f  +  I  tan  9)v  "  IcTsT  a    (13) 


3v       u   9v   v  9v    ,c    .    u  .    ftv     1  W      /n/s 

H=-^^a-7^-  (f +7tan6)u- al?  (u> 


It  =  -  irb  {fr<*u)  +  !i<*v  cose»  <15> 


and  the  flux  form  of  (13)  and  (14) 

^-  -  -  -db^l7<*u2>  +  !«<* uv  cose>}  +  (f  +  *  tane)*v  -  7^lr^2/2) 

3t       acos8  3A        38  a  acos8.  a\ 

(16) 


ft"(*v)  =  "  Ic^?{i(*  UV)  +^(*v2cos6)}  -  (f  +^  tan8)»u  -  \  ^ (*2/2) 

(17) 


In  their  numerical  scheme  they  use  a  centered  time  difference  and  a 
centered  space  difference  for  all  derivatives.   For  the  terms  involving 
(f  +  —  tan8)   they  use  a  time  average   (u  =  —  [u(t+l,  x)  +  u(t-l,  x)]) 
so  that  these  terms  in  either  (13)  and  (14)  or  (16)  and  (17)  are  treated 
implicitly  -  thus  creating  a  semi-implicit  model. 

They  ran  three  tests.   In  the  first,  they  used  5°  grids  in  both  A 
and  8,  taking  the  time  step  based  on  the  smallest  spacial  distance  near 
the  pole.   In  the  second  they  arbitrarily  omitted  points  on  longitude 
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circles  as  the  latitude  increased  towards  the  pole,  so  that  at  60°N  they 
had  72  points  per  longitude  circle  and  at  85°N  they  had  12.   (This  has 
the  same  effect  as  an  icosohedral  grid).   This  allows  a  longer  time  step, 
by  a  factor  of  6.   In  the  third,  they  took  the  output,  did  a  longitudinal 
Fourier  analysis  and  threw  away  the  high  frequency,  short  wave  components 
("filtering").   (This  is  what  is  proposed  for  FNWC  P.E.  model).   They 
start  up  the  scheme  with  a  known  steady  state  solution  of  (13)  -  (15), 
namely  zonal  geostrophic  flow,  given  by 


u  =  u  (cos8  cosa  -  cosA  sin6  sina) 
o 

v  =  u  sinA  sina     0  (18) 

o  2 

Uo  2 

<j>  =  <j>  -  (a  ft  u  H — )  (cosA  cos8  sina  +  sine  cosa) 


(If  a  =  it/2  there  is  no  tendency  to  flow  across  the  equator,  which 
simplifies  the  calculation) . 

Their  results  are  -  for  test  one 

1)  For  the  initial  time  step  the  error  in  u  is  dominated  by  the 
truncation  error  in  $. 

2)  The  error  in  u  grew  by  two  orders  of  magnitude  in  five  days, 
from  the  error  after  one  time  step. 

3)  The  errors  start  at  the  pole  and  spread  out. 
In  test  two,  there  was 

1)   a  factor  of  ten  in  the  size  of  the  error  near  the  poles. 
In  test  three  they  used  the  same  time  step  as  in  test  two  (six  minutes 
vs  one  minute  for  test  one).   The  computations  were  three  times  faster 
than  in  test  one.   They  were  able  to  reproduce  exactly  the  results  of 


test  one.  With  a  fourth  order  difference  scheme  in  space  and  the  filtered 
scheme  they  were  able  to  get  the  best  results  by  about  two  orders  of 
magnitude.   Thus  for  this  type  calculation  a  filtered  model  is  clearly 
the  best. 

Kwizak  and  Roberts  [3]  rewrite  the  three  equations  (13)  -  (15)  as 
follows: 

Let 

K  e  i  (u2  +  v2)  (19) 

and 

la      a 

(20) 


QEF  +  ic^fiv-^ucose>i- 


Then  we  can  write 


_3u 
3t 


acos8  3A 


!£  + 


Qv  - 


3K 


acos9  3A 


(21) 


3v  _  _  1  li  _ 


3t 


a  ;    ^U   a  39 


(22) 


3t 


7^{lx(*u)  +  I?((()V  cos6)} 


(23) 


The  function  K  is  interpretable  as  the  kinetic  energy  and  Q  is 
interpretable  as  the  absolute  vorticity  [3]  or  the  potential  vorticity 
[4].   The  latter  comes  from  the  fact  that  (20)  can  be  written 


Q  =  F  +  V  x  V 

where  V  is  the  velocity  vector. 

Now  the  original  equations  are  derivable  from  assuming  an 
incompressible  gas,  so  that  we  have  the  implicit  equation 


(24) 
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■^  +  ^ =  °  (25) 


which   implies    that   there   is   a  stream  function     ty     such   that 


3d>  3u> 

u  =  -  — ■*-  .      v  =  — x 

3y   *  3x 


(26) 


so  that  (23)  can  be  rewritten 


V24j  =  Q  -  F  (27) 


And  taking  the  curl  of  (21)  -  (22)  we  get 


9t   acosB  3X   a  36  U  ; 


which  can  be  rewritten 


|£-J(Q.llO  (29) 


where  J  is  the  Jacobian  of  Q  and  ip.   The  pair  (26)  and  (29)  are  a 

coupled  system,  with  (26)  being  an  elliptic  equation,  which  could  be 
solved  iteratively,  and  then  <j>  could  be  found  from  (15). 

Observe  that  taking  the  curl  of  (21)  -  (22)  eliminates  the  cj>  from 

these  equations.   Typical  values  of  the  variables  in  the  atmosphere  are 

4 
(j)  ~  3  x  10  ,  u  =  5,  v  =  1.   So  a  potentially  large  term  has  been 

eliminated  from  the  computations  for  the  velocity  V.   Also  observe  that 

—ft 

(28)  is  the  only  wave-like  equation  left  for  V  and  this,  being  first 
order,  has  only  one  wave  solution.   A  linearization  of  (28)  about  a 
steady  flow  V   will  give  only  one  wave 

ik(x-V  t) 


.Q_ 


which  is  the  "meteorological"  wave  which  is  expected:   that  is,  the  two 
gravity  waves  do  not  enter  into  the  computation  for  Q.   This  is  the 
reason  that  the  vorticity  form  of  the  equations  was  used  for  many  years  - 
the  time  dependent  portion  of  the  problem  involved  only  a  single  wave 
speed,  the  one  which  was  desired.   The  problem  is  that  the  solution  must 
be  kept  divergence  free;  this  may  be  a  nontrivial  task.   For  experimental 
purposes  we  can  see  that  the  "meteorological  wave"  ought  to  be  associated 
with  having  a  divergence  free  field,  and  thus  making  an  effort  to  get  the 
initial  conditions  divergence-free  is  a  worthwhile  computation,  as  far 
as  reducing  "unwanted  noise"  in  the  answer.   For  this  reason  Kwizak  and 
Roberts  comment  "...  the  winds  are  perfectly  non-divergent'  initially  and 
at  the  end  of  the  first  time  step.   This  property  virtually  eliminates 
the  gravity  waves  from  the  integration,"  [3]. 

The  computations  which  Kwizak  and  Roberts  do  are  based  on  the  three 
primitive  equations  (21),  (22)  and  (15),  doing  a  semi-implicit  scheme. 
They  take  the  <j>  terms  in  (21)  and  (22),  and  then  rewrite  (15)  as 

and  then  treat  the  first  term  implicitly.   They  thus  have  three  coupled 
equations  to  handle.   The  actual  method  of  computation  is  to  observe 

that  (21)  and  (22)  can  be  used  to  convert  (29)  to  a  non-homogeneous 

-2t  -2t 

Helmholtz  equation  for  <j>  '  ,   defined  as   <f>    e  [<J)(t+At)  +  <j)(t-At)], 


namely 


^(At)2*/*"2'  -  <f 2C  =  f«.».»)  (31) 
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Having  found  this,  then  (21)  and  (22)  can  be  used  to  compute  u  and  v 
in  a  purely  explicit  fashion.   For  their  computations  they  can  use  a  60 
minute  time  step  (vs  a  10  minute  for  an  explicit  calculation)  at  a 
savings  of  4  times  in  the  computation  speed.   Thus  the  method  of 
computation  does  not  involve  three  coupled  implicit  equations  but  one 
Helmholtz  equation  (second  order)  and  two  explicit  calculations. 

Two  years  later  Elviers  and  Sundstrom  [5]  do  a  similar  test,  using 
essentially  the  same  scheme,  noting  that  the  averaging  operator  which  is 
used  allows  for  a  decoupling  of  the  equations  for  even  and  odd  time  steps s 
i.e.,  a  staggered  grid.   They  do  a  stability  analysis  of  the  finite 
difference  models.   Their  analysis  shows  the  semi-implicit  scheme  is 
superior  to  an  explicit  scheme. 

Williamson  and  Browning  [2]  and  Kwizak  and  Roberts  [3]  both  do  a 
semi- implicit  scheme,  but  they  differ  as  to  which  terms  they  treat 
implicitly.   Kwizak  and  Roberts  do  the  more  usual  method.   In  the  u  and 
v  equations  they  treat  the  pressure  gradient  term  (involving  <j>) 
implicitly.   This  allows  them  to  eventually  eliminate  u  and  v  from 
the  <f>  equation  and  convert  the  latter  to  a  second  order  Helmholtz 
equation.   This  one  equation  is  solved  fully  implicitly,  with  u  and  v 
then  found  explicitly.   Williamson  and  Browning  treat  the  Coriolis  terms 
implicitly,  giving  a  coupled  system  for  u  and  v,  while  the  <j> 
equation  is  treated  explicitly.   In  both  cases  only  one  of  the  terms 
driving  the  gravity  waves  is  treated  implicitly,  but  this  is  sufficient 
to  remove  them  from  the  Courant  criterion. 

Thus  a  significant  problem  for  the  predictive  equations  is  the 
ability  to  solve  quickly  and  efficiently  a  system  of  partial  differential 
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equations.   The  usual  method  is  to  use  classical  finite  difference 
techniques.   The  standard  way  is  to  use  second  order  centered  differences. 
There  is  every  reason  to  believe  that  fourth  order  centered  differences  in 
space  will  increase  the  accuracy.   This  is  currently  being  investigated  [6] 

3.   Finite  Element  Methods 

A  question  of  definite  interest  is  whether  the  "finite  element 
method"  of  solving  partial  differential  equations  will  give  better 
accuracy.   I  am  aware  of  only  four  places  where  this  is  currently  being 
investigated.   George  Fix  [4,7]  is  studying  ocean  circulation  problems 
this  way.   His  studies  are  being  continued  by  Hirsch  [17].   M.  P.  Cullen 
[8,  9,  13]  has  programmed  the  barotropic  equations,  and  is  now  attempting 
to  program  a  more  realistic  set  of  primitive  equations.   A.  Staniforth 
[10]  is  attempting  to  implement  finite  element  calculations  in  the 
Canadian  Meteorological  Office.   And  a  student  thesis  at  NPS  by  Donald 
Hinsman  [11]  ran  some  experiments  with  finite  elements  on  the  barotropic 
equations.   All  of  these  indicate  that  this  method  may  have  a  significant 
future  in  meteorology. 

Fix  [4]  looks  at  the  ocean  circulation  problem,  which  is  just  (13) 
and  (14) ,  together  with  the  divergence  free  condition 

V  •  V  =  0 

Thus  he  does  not  get  involved  with  gravity  waves  and  has  only  the 
"meteorological  mode"  to  contend  with.   He  then  converts  to  the  vorticity 
equation  and  the  non-linear  balance  equation  (27),  (29).   He  then  takes  a 
finite  element  -  Galerkin  approach,  using  linear  elements  (also 
quadratics  and  cubics  for  further  tests). 
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There  are  three  problems  to  be  addressed  in  any  analysis  of  this 
discretization.   The  first  is  the  accuracy  of  the  spacial  discretization. 
The  second,  which  is  due  to  the  fact  that  these  equations  are  non-linear, 
is  called  "aliasing,"  a  feature  which  can  be  most  easily  seen  by 
considering  a  Fourier  analysis  of  the  spacial  terms.   If  both  u  and  v 
are  written  as  Fourier  series  and  then  truncated  at  the  N —  harmonic, 
then  a  term  involving  u  times  v  will  have  a  term,  say 
(u^cos  N  x) (v  cos  M  x),  which  would  normally  give  rise  to  two  terms 
2  ^Vm  cos(m-N)x  anc*  ~y   ^"^vr  cos(N+M)x,   but  the  latter  can  not  appear 
due  to  the  truncation.   Certain  discretization  schemes  have  an  imbalance 
in  the  treatment  of  this  phenomenon,  known  as  aliasing.   However, 
Jesperson  [12]  has  shown  that  this  phenomenon  does  not  occur  with  a 
finite  element  scheme,  that  is,  finite  element  schemes  are  free  from 
aliasing,  a  fact  which  Fix  reconfirms. 

Fix  also  shows  what  is  widely  known,  that  the  spacial  accuracy  for 
the  velocity  is  0(h  )  where  h  is  the  grid  size  and  k  is  1,  2,  or 
three  depending  upon  whether  linear,  quadratic  or  cubic  elements  are 
used. 

The  third  problem  is  how  to  handle  the  time  integration.   Fix  does 
not  have  to  worry  about  a  semi-implicit  scheme  from  the  point  of  view  of 
gravity  waves.   However,  any  finite  element  scheme  links  more  than  two 
points  and  one  is  automatically  forced  into  an  implicit  scheme.   That  is 
one  drawback  for  finite  elements.   Fix  proceeds  to  analyze  a  linear  one 
dimensional  analogue  of  the  wave  propogation  problem.   The  finite 
element  method  generates  its  own  natural  set  of  difference  equations. 
For  the  linear  problem  which  Fix  sets  up,  using  linear  elements  in  space 
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with  time  varying  coefficients,  Fix  shows  that  the  implicit  discretization 
which  the  finite  element  method  forces  gives  fourth  order  accuracy  in  the 
phase  speed  of  waves,  a  fact  which  Cullen  had  noted  earlier  [8]. 

Fix  then  chooses  for  his  time  discretization  the  usual  centered 
difference  (leap-frog)  time  discretization,  so  that  the  final  computation 
for  Q  is  given  by 

jjf[Qn+1   -  (f-^.dA  =  2At    Jjfj^'    *"><*  • 
A  A 

Thus  for  the  non-linear  ocean  circulation  problem  one  knows  that 
the  finite  element  method  does  not  introduce  aliasing,  is  spatially  as 
accurate  as  we  make  the  finite  elements,  and  can  be  conjectured  to  be 
fourth  order  accurate  in  phase  speed.   (Phase  speed  of  the  waves  has 
always  been  a  problem  in  forecasting  of  weather). 

In  a  series  of  three  papers  [8],  [9],  [13],  Cullen  tackles  two 
problems.   The  first  problem  is  to  solve  equations  (3)  in  a  limited  area 
-1  £  x,  y  <_  1   (the  "beta  plane")  with  periodic  boundary  conditions, 
using  a  finite  element  method.   The  second  is  to  solve  equations  (3)  on 
the  surface  of  the  globe.   His  analysis  proceeds  as  follows. 

He  first  [13]  considers  a  single  linear  equation 

||  +  V-VcJ)  =  0  (32) 

on  a  rectangle,  using  rectangles  and  bilinear  elements,  with  V  known. 
He  compares  his  results,  using  a  16  x  16  grid,  to  second  and  fourth 
order  finite  difference  schemes  using  a  32  x  32  grid.   The  exact  solution 
to  this  problem  is  known.   He  demonstrates  that  the  finite  element 
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calculation  is  more  accurate,  both  in  handling  smooth  data  and  in  handling 
a  problem  with  a  discontinuity. 

Satisfied  with  the  results  of  the  linear  problem,  Cullen  [13] 
proceeds  to  the  non-linear  problem  (3),  in  a  grid  -1  <_  x,  y  <_  1,  without 
the  Coriolis  term.   His  initial  condition  would  give  a  gravity  wave  if 
the  problem  were  linear.   He  compares  his  answer  with  a  16  x  16  grid  to 
finite  difference  schemes  with  32  x  32  and  64  x  64  grids.   He  uses  linear 
finite  elements  with  a  leap-frog  time  step.   His  results  indicated  that 
the  finite  element  method  was  better,  although  there  is  no  exact  solution 
to  compare  with.   Cullen  then  attempts  to  analyze  the  numerical  scheme. 
For  the  linear  one  dimensional  case  he  shows  that  the  phase  error  is 
fourth  order  in  At/Ax.   In  fact,  he  shows  something  more,  namely  that 
for  some  range  of  At/Ax  the  phase  error  can  actually  cause  a  small 
leading  phase.   He  also  argues,  largely  on  qualitative  grounds,  that 
there  is  no  aliasing  problem.   There  is  no  discussion  of  the  inertial  vs 
Rossby  waves,  and  no  attempt  to  isolate  one  from  the  other,  or  to  control 
either. 

Cullen' s  second  paper  [8]  concerns  equations  (3)  on  the  entire 
globe  -  the  genuine  barotropic  problem.   One  question  which  Cullen  now 
addresses  is  how  is  the  best  way  to  handle  the  non-linear  terms  which 

appear  in  (3).   His  approach  is  to  analyze  a  single  one  dimensional  term 

9v 
of  the  form  u  —  .   The  various  possibilities  are: 
3x 

1.  Treat  v  as  known,  take  the  nodel  values  of  v,  compute  a 
derivative,  and  simply  multiply  the  coefficients  for  u  by  these  values. 

2.  Compute  a  finite  element  expansion  for  v,  analytically 

8v 
differentiate  this  expansion  to  find  —  ,  and  use  the  resultant  two 

oX 

expansions  to  compute  t  =  u  —  . 
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9v 
3.   Compute  a  finite  element  expansion  for  —  itself. 

Cullen  claims  that  the  last  method,  while  sacrificing  some  accuracy, 

controls  aliasing  much  better  than  either  of  the  first  two. 

I  was  unable  to  follow  his  arguments,  and  will  attempt  in  a  further 
report  to  see  whether  it  can  be  generalized.   In  any  event  Cullen  uses 
the  third  method  in  his  calculations  [8]. 

In  [8]  Cullen  runs  four  tests  on  the  equations  (3),  using  the  three 
methods  above  (for  two  of  the  runs  he  slightly  modifies  the  coefficients 
in  method  3) .   He  takes  as  initial  values  a  finite  element  solution  to 
(2a),  comparing  with  other  published  results  using  finite  differences. 
Assuming  that  the  published  finite  difference  results  are  the  most 
accurate  (?)  he  notes  that  the  phases  on  most  finite  element  runs  appear 
to  lag,  although  some  are  advanced.   He  concludes  that: 

(a)  Waves  down  to  four  element  lengths  will  be  treated  almost 
perfectly. 

(b)  Waves  less  than  two  grid  lengths  will  not  be  treated  at  all  well, 

(c)  The  finite  element  scheme  "essentially  eliminates  aliasing 
errors." 

(d)  The  boundary  condtions  used  on  the  problem  introduced  errors  of 
the  same  order  of  magnitude  as  the  change  from  finite  differences  to 
finite  elements,  and  are  thus  quite  significant. 

In  his  third  paper  [9]  Cullen  reports  on  actual  computations  on  a 
sphere,  relying  heavily  on  the  analysis  above.   He  concludes  "the  finite 
element  method  is  computationally  more  efficient."  He  uses  an  icosahedral 
grid,  subdivided  by  latitude  and  longitude  lines  to  form  triangles, 
resulting  in  1002  points.   In  integration  he  treats  the  trig  functions  in 


-16- 


(3)  as  constants  over  each  small  triangle.   He  uses  the  scheme  of  [8]  to 
treat  the  non-linear  terms. 

For  his  actual  computational  scheme  it  appears  that  the  only  implicit 
portion  of  the  scheme  is  the  implicitness  generated  by  the  left  hand  side. 
In  other  words,  looking  at  (3)  as  equations  of  the  form 

he  solves  the  resultant  non-linear  equations  by  the  following  iterative 
process.   Treat  the  right  hand  side  as  known  and  take  the  diagonally 
dominant  matrix  on  the  left  as  the  generator  of  the  next  iteration.   In 
this  way  the  iterative  process  to  find  the  values  of  the  nodal  points  at 
time  t  is  relatively  fast.   Then  a  leap-frog  time  discretization  of  10 
minutes  was  used.   (He  could  use  one  of  up  to  14  minutes).   He  found  that 
filtering  was  required  every  two  hours  to  get  long  time  (greater  than  5 
days)  solutions.   His  initial  conditions  were  Rossby  waves  with  wave 
number  4  and  wave  number  8.   He  compares  his  results  with  published 
results  using  a  finite  difference  model  with  4032  points,  one  with  14,592 
points  and  a  spectral  model  with  640  degrees  of  freedom.   His  results  are 
better  than  the  first  but  not  as  good  as  the  last  two.   He  also  observes 
that  the  errors  seem  to  start  from  the  vertices  of  the  icosahedron,  where 
there  are  only  5  supporting  triangles  instead  of  6  as  there  are  at  the 
intermediate  triangulation  points. 

Hinsman  [11]  in  his  master's  thesis  again  studied  equations  (3).   He 
considered  two  possible  grids.   One  used  the  two  angles  X  and  9   as 
rectangular  coordinates  and  triangulating  the  resultant  "rectangle." 
This  results  in  a  very  fine  subdivision  of  the  polar  regions  resulting  in 
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very  small  physical  lengths  Ax  in  that  region.   The  second  grid  was  an 
icosahedral  grid,  (as  did  Cullen  [8],  [9])  which  was  subsequently  sub- 
divided by  arcs  of  great  circles.   This  results  in  most  of  the  nodes 
having  6  adjacent  nodes.   The  corners  of  the  icosahedron  have  only  5 
adjacent  nodes.   This  appears  to  generate  some  "noise,"  as  Cullen  noted 

[9]. 

Experiments  with  the  (X,6)  grid  showed  instability  after  12  hours, 
as  might  be  guessed  from  the  fact  that  there  were  36  points  around  each 
latitude  circle  including  those  of  the  poles.   The  instability  clearly 
arose  in  the  polar  region.   Experiments  with  the  icosahedral  grid  did  not 
show  these  problems. 

Starting  with  an  analytic  solution  to  the  non-linear  balance  equation 
which  has  essentially  one  wave  going  around  the  earth,  Hinsman  Fourier 
analyzed  the  solution  as  it  propagated,  comparing  his  results  to  similar 
finite  difference  calculations.   In  the  low  latitudes  the  propagation 
speed  was  almost  exactly  correct;  as  compared  to  50-60%  correct  for 
finite  differences.   In  the  high  latitudes  both  finite  elements  and 
finite  differences  fall  behind  the  predicted  speed.   (This  contradicts 
Cullen' s  observations). 

The  method  of  solution  was  very  different  from  any  previous  methods. 
Each  equation  in  (3)  is  "quasi-linearized"  by  considering  the  other  two 
variables  as  known  (from  a  previous  calculation).   The  finite  element 
scheme  automatically  generates  an  implicit  scheme,  but  the  decision  was 
made  to  go  one  step  further  and  use  a  Crank-Nicholson  approach, 
essentially  calculating  the  variables  at  time  step  (N+l/2).   The  quasi- 
linearization  uncoupled  the  equations,  but,  like  Cullen,  there  was  no 
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attempt  to  distinguish  between  the  Rossby  wave  and  the  gravity  waves. 
Finite  elements  were  used  to  expand  all  the  variables,  having  first 
integrated  by  parts  to  get  the  weak  form.   The  equations  were  then 
successively  solved,  first  for  <j>,  then  for  u  and  finally  for  v. 
(The  trig  functions   sin8   and  cos 8  were  also  written  as  finite 
elements,  as  opposed  to  Cullen  who  treated  them  as  constants  over 
triangles).   The  resultant  equations  were  solved  by  a  Gauss-Seidel 
iteration.   This  took  about  10-12  iterations  to  converge.   (An  additional 
feature  of  the  program  was  a  very  efficient  coding  scheme  to  avoid 
completely  the  storage  of  the  zero  elements  of  the  matrices  involved). 

Swartz  and  Wendroff  [14]  compare  the  relative  efficiency  of  finite 
difference  and  finite  element  methods.   They  do  so  by  taking  a  one 
dimensional  linear  problem.   Their  interest  is  somewhat  different  than 
ours  and  thus  their  data  is  recorded  very  differently.   They  record,  for 
example,  (Table  1)  the  number  of  intervals  per  wave  length  which  are 
necessary  to  get  a  desired  accuracy  in  phase  (in  theory).   One  conclusion 
which  can  be  drawn  is  that  a  dramatic  improvement  in  resolution  can  be 

obtained  by  switching  from  linear  elements  to  quadratic  elements  (for  an 

-2 
error  of  10  '  the  number  of  intervals  per  wavelength  goes  from  8.7  to 

-4 
4.8,  while  for  10   error  the  change  is  from  27  to  9.7). 

They  also,  using  results  of  Kreiss  and  Oliger,  report  the  results  of 

finite  difference  theory.   For  fourth  order  spacial  accuracy,  -  which  is 

what  linear  elements  give  -  the  corresponding  results  are  a  change  from 

13.3  to  7.9  and  from  42.5  to  17.3.   So  one  can  conclude  that  in  theory 

the  finite  elements  have  better  phase  resolution  (e.g.,  8.7  intervals  per 

wave  length  vs  13.3  for  linear  elements  and  4th  order  differencing)  and 

that  the  payoff  for  increasee  complexity  is  also  greater. 


-19- 


The  method  that  they  use  for  actually  computing  the  finite  element 
solutions  is  interesting  -  they  call  it  "use  of  the  trapezoidal  rule." 
Consider  the  equation 


du  3u 

3t     °  9x 


and  compute  the  finite  spacial  element  expansion  of  this  to  get  a  formula 

A  u  =  B  u  (33) 

where  A  and  B  are  sparse  matrices  and  u  is  a  vector  (for  linear 


elements  A  has  a  —  (1,  4,  1)  structure).   Now  take  a  time  discretization 


which  gives 


~n+l   ~n      ~n+l  .  ~n 
A^^  -  B^^  (34) 


This  scheme  is  fully  implicit  and  is  stable  for  all  h,   At.   They  compare 
this  with  the  leap-frog  finite  difference  scheme. 

Before  finishing  with  the  linear  problem  they  conclude  that  a 

twelfth  order  spacial  difference  scheme,  with  leapfrog,  is  competitive 

-4 
with  a  cubic  spline,  for  phase  errors  of  10 

They  turn  to  the  first  order  one- dimensional  non-linear  problem 

u  -  |-  f (u)  +  g(u)  .  (35) 

They  have  two  conclusions,  one  of  which  is  given  just  in  passing.   They 

remark  that  they  will  not  use  finite  elements  directly  on  (41).   Instead 

-1  8 

they  will  use  the  operator  A  B  to  replace  —  ,   where  A  and  B   are 

oX 

defined  in  (33).   Their  reasoning  is  that  the  finite  element  scheme  is 
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awkward  and  that  it  has  been  shown  that  use  of  the  finite  element  scheme 
directly  will  degrade  the  spacial  accuracy  [15].   The  scheme  above  will 
keep  the  truncation  error.   They  further  remark  that  "It  no  longer  makes 
sense  to  try  to  find  the  best  scheme.   A  reasonable  approach  seems  to  be 
to  compare  schemes  using  the  same  number  of  intervals  per  wave  length." 
This  sort  of  reasoning  is  directly  counter  to  what  Cullen  and  Hinsman 
have  used. 

Their  second  conclusion,  based  on  the  above  and  counting  evaluations 
of  f(u)   and  g(u)   as  the  major  source  of  computer  time.,  is  that,  using 
cubic  splines  and  18th  order  differencing  as  comparable  schemes,  finite 
difference  and  finite  element  schemes  are  roughly  comparable.   (For 
global  meteorology  an  18th  order  difference  scheme  seems  excessive) . 

Another  paper  which  deals  with  non-linear  finite  element  calculations 

for  hyperbolic  problems  is  by  Oden  and  Fost.   What  Oden  and  Fost  do  to 

find  the  finite  element  formula  for  a  hyperbolic  problem  is  to  use  as  a 

test  function   V.  I  t—  J  ,   so  that  the  equation  with  which  they  work  is 
1'  VI  at 


)' 


dE 
the  analogue  of  -7-  =  0  ,  where  E  is  the  energy  of  the  system.   It  is 

then  clear  that  their  method  preserves  energy.   In  the  actual  computation 

3V. 
they  use  the  centered  (leapfrog)  approximation  to  —7—  ,  and  their 

o  t 

analysis  is  of  a  single  second  order  equation 

2  2 

3  u     3T  t      .  3  u   n  f„. 

-   c  - — (u  )  — r  =  0  ,  (36) 


„  2     3u  v  x'  .  2 

3t        x     3x 


which,  in  weak  form,  is 


32u   3u\     2  (l,      s      32u 
^2  '  -oil  +  C  [T(\^   ix^ 
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Their  study  reveals  the  following: 

h    fl 

(a)  the  resultant  equation  is  numerically  stable  if  —  >^—  C,   N , 
v  '  n  J  At    3   (max) 

where  C,        s      is  the  maximum  wave  speed  C^T' (u  )  . 

(max)  x      .  . 

(b)  If  instead  of  using  the  matrix  generated  by  | — r  ,  —  |  (the 

V9t     7 

so-called  consistent  mass  matrix)  one  uses  a  diagonal  matrix  (the  lumped 
mass  approximation)  the  stability  criterion  is  improved  to  read 

At       (max) 

(c)  The  finite  element  solution  converges  uniformly  to  the  solution, 

using  piecewise  linear  elements.   They  do  not  make  any  analysis  of  the 
order  of  accuracy  of  this  approach,  so  this  paper  does  not  touch  the 
question  raised  by  Swartz  and  Wendroff. 

4.   Questions  and  Some  Answers 

A  number  of  questions  arise  from  the  intersection  of  these  papers. 

1.  What  is  the  best  method  of  handling  the  non-linear  terms,  and 
the  variable  coefficients?   Cullen,  Hinsman,  Oden  and  Swartz  all  advocate 
different  answers  for  different  reasons. 

2.  What  does  the  finite  element  method  do  to  the  gravity  waves 
which  worry  the  people  doing  finite  difference  calculations?   None  of  the 
finite  element  calculations  even  mention  this  question. 

3.  What  are  the  merits  of  the  "lumped  mass"  vs  "consistent  mass" 
approach  which  is  so  familiar  to  the  mechanical  engineers? 

4.  If  this  were  to  be  implemented  as  an  operational  scheme,  what 
are  the  possible  merits  of  SOR  or  ADI.   Hirsch  [17],  doing  ocean 
circulation  problems,  solves  the  Poisson  equation  (27)  by  SOR  and  the 
advection  equation  (29)  by  ADI.   Staniforth  and  Mitchell  [10]  use  an  ADI 
method  for  their  calculations. 
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5.  Cullen  does  not  integrate  by  parts  to  compute  his  u  — 

oX 

term.   What  effect  does  this  have  on  the  solution?   Hinsman  (unpublished) 
noted  a  considerable  improvement  in  his  results  when  he  used  the  weak 
form  (integrated  by  parts)  of  the   $   equation  as  opposed  to  the  strong 
form  (3) . 

6.  Fix  has  proved  that  the  vorticity  equations  for  ocean  circulation 
automatically  satisfy  the  desired  conservation  laws  when  written  as 
finite  elements.   Is  this  also  true  for  the  barotropic  equations? 

7.  What  does  the  finite  element  method  do  to  the  phase  speed  in 
two  dimensions? 

We  now  proceed  to  answer  this  last  question;  at  least  for  a  linear 
model.   The  finite  element  approximation  to 

|J+  V-V<j>  =  0  (37) 

using  bilinear  elements   on  rectangles   is 


f {16   -h.  +  4[iM+1  +  lKi_x  +  im>.  +  i,_w]  +   [*m>.+1  +  Vlij+1  + 

♦m-1,1-1  +  Vi,j-i1}  +  36s"4  W.j  -  4  Vi.j1  +  [*m,i+i '  Vi.j+i1 

We  investigate  the  phase  propagation  of  this  linear  model.   The  exponen- 
tial solution  to  (37)  is  $  =  &   e-i(wt-  *x)  where  u  =  y.k>   If  we  look 

* 

for  u   in  the  approximation 
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* 


then 


thus 


a  .*  =  A  e_iw  t  +  Uk+ijh  (38) 


-ikhw  {16  +  8 [cos  h  +  cos  k]  +  4  cos  k  cos  h} 


+  3hui{sin  k}(2  +  cos  h)  +  3kvi{sin  h}(2  +  cos  k)  =  0 


*    3uh  sin  k(2  +  cos  h)  +  3vk  sin  h(2  +  cos  k) 
w  hk(2  +  cos  h)(2  +  cos  k) 


3  sin  k  3  sin  h  ,  ~qv 

u  (2  +  cos  k)k       (2  +  cos  h)h  Uy; 


We  then  have  the  following  theorem. 

Theorem  1:   The  finite  element  approximation  to  (37)  using  bilinear 

elements  on  rectangles  has  fourth  order  accuracy  in  phase  speed. 

Proof:      From   (39)  ?  , 

k  k 

(3-|  +3i  ~  •••  > 

M  .  u  —  —  

(3-|  +  i     -    ...  ) 

h2      >,4 

(3  -f    +  3  I       -      ...    ) 

+     


h2         h4 
^       2  4!  ; 


/.  w*  ~   u(l  +  0(k4))  +  v(l  +  0(h4))  (40) 

Q.E.D. 


This  is  to  be  expected,  since  the  basis  elements  are  the  tensor  product 
of  two  one- dimensional  linear  elements  on  rectangles. 
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A  more  interesting  result  is  whether  the  same  phenomenon  happens  if 
we  use  linear  elements  on  triangles.   We  have  the  following  result.   The 
finite  element  approximation  for  (37)  using  the  general  linear  element 
<J>  =  a  +  bx  +  cy  on  triangles  is 


kh   •      •        *        •        •        •  • 

T2{6  h,i  +  Vi.j  +  ♦t.j+i  +  Wj  +  ♦t.j-i  +  ♦m.j+i  +  ♦i-i,j-i}  + 

vk 


This  also  has  a  fourth  order  phase  speed  accuracy  as  we  now  show.   If  we 
again  look  at  $  *  as  in  (38)  we  get,  for  the  first  term  in  (41) 

* 

-  — rr —  {3  +  cos  k  +  cos  h  +  cos  k  cos  h  -  sin  k  sin  h} 

2ikhco*   t,       .2       ,2       . ,     ,    hk3   ,    kh3   ,  , 

+** ttj —  {6-k     -h     -kh+-rr     +  Tj~     +•••    J 

while  the  second  term  is 

-^7^  {3  sin  k  -  sin  h(l  -  cos  k)  -  sin  k(l  -  cos  h) } 
6 

,  2ikhu  ,,   .2   ,2   ..  ,  h_k       , 
r«J  — — —  {6-k  -h  -  kh  H 7-  +  . . .  > 

and  the  third  term  is 

^~  {3  sin  h  -  sin  k(l  -  cos  h)  -  sin  h(l  -  cos  k) } 
6 

_.  2ikhv  r,       ,2   ,2       k  h       -, 
r*»  — — —  {6-k  -h  -  kh  +  —  +  . . .  t 
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Thus 

w*~  u(l  +  0(h\n))  +  v(l  +  O(hV))  (42) 

where  m  +  n  =  4.   Thus  we  have  shown 

Theorem  2:   The  finite  element  approximation  to  (37)  using  linear 
elements  on  triangles  has  fourth  order  accuracy  in  phase  speed. 
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